Dynamical response of a neuron–astrocyte coupling system under electromagnetic induction and external stimulation
Yuan Zhi-Xuan, Feng Pei-Hua, Du Meng-Meng, Wu Ying
State Key Laboratory for Strength and Vibration of Mechanical Structures, Shaanxi Engineering Laboratory for Vibration Control of Aerospace Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, China

 

† Corresponding author. E-mail: wying36@xjtu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 11772242) and China Postdoctoral Science Foundation (Grant No. 2018M631140).

Abstract

Previous studies have observed that electromagnetic induction can seriously affect the electrophysiological activity of the nervous system. Considering the role of astrocytes in regulating neural firing, we studied a simple neuron–astrocyte coupled system under electromagnetic induction in response to different types of external stimulation. Both the duration and intensity of the external stimulus can induce different modes of electrical activity in this system, and thus the neuronal firing patterns can be subtly controlled. When the external stimulation ceases, the neuron will continue to fire for a long time and then reset to its resting state. In this study, “delay” is defined as the delayed time from the firing state to the resting state, and it is highly sensitive to changes in the duration or intensity of the external stimulus. Meanwhile, the self-similarity embodied in the aforementioned sensitivity can be quantified by fractal dimension. Moreover, a hysteresis loop of calcium activity in the astrocyte is observed in the specific interval of the external stimulus when the stimulus duration is extended to infinity, since astrocytic calcium or neuron electrical activity in the resting state or during periodic oscillation depends on the initial state. Finally, the regulating effect of electromagnetic induction in this system is considered. It is clarified that the occurrence of “delay” depends purely on the existence of electromagnetic induction. This model can reveal the dynamic characteristics of the neuron–astrocyte coupling system with magnetic induction under external stimulation. These results can provide some insights into the effects of electromagnetic induction and stimulation on neuronal activity.

1. Introduction

Neurons are well known as the main elements of the central nervous system (CNS). The Hodgkin–Huxley neuron model[1] opened up the research of neurons from the viewpoint of dynamics. Although both the Morris–Lecar model[2] and FitzHugh–Nagumo (FHN) model[3] are two-dimensional simplifications of the Hodgkin–Huxley model, they can also correctly represent the firing modes of neurons. The Hindmarsh–Rose model,[4] another simplified Hodgkin–Huxley model, can be employed to show neuronal discharging properties near the threshold potential without external stimulation, such as noise, instead of strong stimulation. The firing properties of neurons can also be analyzed by bifurcation theory.[58] Some neuron models[912] have been established for further research in computational neuroscience. Gu[13] found that the inhibitory coupling current of delay modulation can cause multiple synchronizations. Ma[14] introduced a new electric field variable into the simplified FHN neuron model, taking into account the effect of the electric field of the neuron. Wang[15,16] observed that the activities of neurons, based on the Hodgkin–Huxley neuron model, have duality that energy expenditure exists in subthreshold neurons and suprathreshold neurons, but only suprathreshold neurons have energy absorption. The duality of neuronal activity has also been found in Chay model[17] and structural neural networks.[18,19] And it has been proved in functional networks such as intellectual exploration.[20,21]

Many physiological experimental results have verified that astrocytes can regulate the transmission of electrical signals among neurons.[2224] Postnov[25] established a detailed model of a tripartite synapse coupling of P-neuron and R-neuron together with a giant astrocyte in the ganglia of the medical leech. Li and Rinzel[26] analyzed and reduced the nine-variable De Young–Keizer model for Ca2+ oscillations in astrocytic calcium stores to a two-variable system, called Li–Rinzel model, and then focused on and revealed the channel gating variables about Ca2+ activation and Ca2+ inactivation. Nadkarni[27] put forward a model for neurons coupled with astrocytes and predicted spontaneous oscillations of seizure-like firings between the neuron and astrocyte without stimulation. Li[28] studied the inhibitory effect of calcium channel blockade in astrocytes on neuronal epileptic firing with a modified GABAergic astrocyte model. Erkan[29] discovered that the astrocyte has a great influence on the neuronal weak signal detection performance, which reveals the stochastic resonance phenomenon relying on the intensity of noise, the detection performance of the neuron increases significantly with the increase of the optimal coupling strength of the astrocyte. Kanakov[30] showed that the astrocytic regulation on neuron may be vital about the information attribute in neuron–astrocyte ensembles. In summary, researches of the nervous system no longer ignore the regulation of astrocytes in brain function or consider the “solo theory” of neurons alone.

Electromagnetic field in the nervous system, induced by the oscillations of neuronal membrane potentials, has been observed in the biological experiment,[31,32] then, an induction current generated from the neuronal electrical activity can feedback to the neuron itself. At present, researchers have proposed that adding the effect of electromagnetic induction to the firing activity of nervous system can be realized by memristor,[3335] which is proposed by Chua first as one of four major electronic components.[36] Wang[18,3739] introduced the effect of electromagnetic induction into the biophysical neuron model by considering the magnetic field produced by neural action potential, which causes violent ion exchange, and demonstrated that the neuronal coding induced by brain activities can be characterised by the theory of energy coding. Lv[34] found that electromagnetic induction spontaneously produced by neuronal alternating depolarization and hyperpolarization is contributive to the memory effect of the nervous system. Feng[40] determined the membrane potential of a single neuron model considered magnetic flux occurs among periodic, quasi-periodic, and chaotic motions, where sharp switching between periodic and quasi-periodic motions appeared in a particular parameter interval. Liu[41] observed the amplitude and frequency of electromagnetic inductions influence the responses of the electrical activities of the hybrid neuron which consists the Hindmarsh–Rose (HR) model with the Wilson model. Wu[42] studied the dynamical responses of HR model in consideration of the effect of electromagnetic induction and discovered the induction current can enhance neuronal bursting activities in contrary to traditional view that the reduction effect provided by the induction should inhibit neuronal electric activities. These aforementioned electromagnetic inductions around neurons are thought to merely affect the feedback of magnetic induction to individual neurons.

There are many studies to investigate the effects of tetanic stimulation on the firing characteristics of neurons. Muramatsu[43] demonstrated that different excitation patterns can be evoked by electrical stimulation on the stimulated cortical layer and acoustic responses are most likely simulated by the electric stimuli in layer 4. Katta[44] found that farther channels activated by stronger stimuli is the main factor responsible for touch sensitivity, rather than higher channel open probability. Liu[41] discovered that the multiple fire patterns of neurons based on a hybrid neuronal model arise and transit successively with the increase of stimulus intensity. Therefore, in addition to the above factors, the influence of tetanic stimulation on neural system is worth studying.

The effect of electromagnetic induction on the neuron–astrocyte coupling system is unknown, and present neuron models rarely address it. Fortunately, a relevant model has been proposed in Ref. [34]. In this paper, we study the dynamic properties of an improved neuron–astrocyte coupled model that takes into account the effect of electromagnetic induction. Our results reveal that different stimulus durations and intensities, including limited or infinite stimulus duration and intensity, can induce different kinds of neuronal electrical activity.

2. Model

The neuron is described by the improved Hodgkin–Huxley model:

where V is the membrane potential and n4 denotes the probability that a potassium ion (K+) channel is open. The probability of opening the sodium ion (Na+) activation gate is m3, and the probability of opening sodium inactivation gate is h. VK, VNa, and VL are the reversal potentials of sodium, potassium, and the leakage system respectively, Cm is the membrane capacitance, and the maximal conductance of the potassium, sodium, and leakage channels are gK, gNa, and gL, respectively. The coefficients in Eq. (2) are given by[28]

Among the other terms in Eq. (1), Iext represents an external forcing current, which will produce multiple dynamic responses to the coupling system with different stimulation durations and intensities, and Imag defines the self-induced electromagnetic current from neuronal activity, given by[34]

where the term k1ρ (ϕ)V describes the negative feedback of magnetic induction on membrane potential when the feedback strength term k1 is positive, while Iastro denotes the additive current generated by the astrocyte, which releases specified amounts of neurotransmitters into the synaptic cleft when the neuron fires. Major neurotransmitters bind to the postsynaptic membrane; however, some can bind to the astrocyte, which can cause the release of IP3 within the astrocyte. The modulation of intracellular IP3 in the astrocyte can be modeled by

where [IP3]* is the concentration of IP3 in the equilibrium state and rIP3 determines the production efficiency of IP3. These values were determined in Ref. [45]. Because of the properties of the Heaviside function (Θ(x)), when the neuron membrane potential exceeds +50 mV, the second production term on the right side will be implemented.

The production of IP3 promotes the release of calcium ions (Ca2+), primarily released from the endoplasmic reticulum (ER). Here, we adopt the Li–Rinzel model[27] to describe the dynamics of [Ca2+]

where q is the activated fraction of IP3 receptors (IP3Rs) on the surface of the ER, Jchannel represents the [Ca2+] that moves from the ER to the inner space of the astrocyte through the IP3R channel, Jpump denotes the [Ca2+] pumped from the inner space of the astrocyte into the ER, and Jleak defines the [Ca2+] that leaks from the ER to the inner space of the astrocyte. These variables from Eqs. (6) and (7) can be expressed as

The coupling of an astrocyte brings about an electric current[28]

where the Heaviside function Θ(x) denotes the astrocytic feedback current to neuron can be activated when the astrocytic calcium concentration exceeds +196.69 nM.

The values of these parameters and constants, as well as their units, are given in Table 1.

Table 1.

Parameters and constants in the simulation.

.
3. Numerical results

In this section, the fourth-order Runge–Kutta algorithm is used to solve the ordinary differential equations above with a time step of 0.01 ms.

3.1. The sensitive dependence of “delay” on stimulus duration

We first focus on the changing of neuronal electric activity in the coupling system with varying the duration of the external stimulus. Several neuronal firing patterns are plotted in Fig. 1, the stimulus intensity is fixed at Iext = 11 μA/cm2. The stimulus durations are (10 s–60 s) for Figs. 1(a)1(f), respectively. In Fig. 1, the first arrow indicates the instant when the external stimulus stops, and the interval from 0 s to this instant is called the stimulus duration (ts); then the second arrow, if any, indicates the instant when the firing of the neuron ceases, and the duration between these two arrows is called “delay (td)”. A magnetic field appears when a neuron generates an action potential. Nevertheless, when a neuron no longer fires and recovers to its resting state, the magnetic field will not disappear right away because of electromagnetic induction.[33,46]

Fig. 1. The firing patterns of the neuron under Iext = 11μA/cm2 with different stimulus durations (ts). The first arrow denotes the instant that the stimulus stops and the stimulus duration (ts) is defined as the interval from 0 s to this instant, while the second arrow, if any, indicates the instant that the neuron recovers to the resting state. The duration between these two arrows is called “delay (td)”. The initial values of neuron membrane potential, astrocytic intracellular calcium concentration, and magnetic flux are 0 mV, 20 nM, and 0 nWb, respectively.

Neuronal firing patterns when ts is approximately 20 s or 50 s are shown in Fig. 2. The values of ts are 19.99 s, 20 s, 20.01 s, 49.99 s, 50 s, and 50.01 s for Figs. 2(a)2(f), respectively; in other words the difference is only 0.01 s from panel (b) to panels (a) and (c) and from panel (e) to panels (d) and (f). However, amazingly, td appears from 0 s in panels (b) and (e) and td = 5.1 s in all four graphs containing td.

Fig. 2. The firing patterns of the neuron under Iext = 11 μA/cm2 with different stimulus durations (ts). ts increases by 0.01 s one by one in both rows and the td is appearing and disappearing alternatively. The initial values are the same as those in Fig. 1.

The phase portrait of the systems with ts = 20 s and ts = 20.01 s, the time series of neuronal membrane potential as shown in Figs. 2(b) and 2(c), are plotted in Fig. 3 for clarifying why the “delay” changes so sharply with a small change in ts. In Fig. 3, the difference between ts in panels (a) and (b) are only 0.01 seconds, but their phase trajectories are quite different. The firing pattern of the neuron in panel (a) enters the resting state directly through the red trajectory after the stimulation stops, while “delay” appears after the stimulation stops in the system in panel (b), which enters the resting state through the blue trajectory after about 5.1 s. The difference of 0.01 s results in “delay” indicating that the orbits into the resting state are selective and sensitive while a very subtle difference in ts can lead the trajectory of the system to enter different orbits.

Fig. 3. The phase diagram of the coupling systems under (a) ts = 20 s, (b) ts = 20.01 s. (a) At the beginning, the neuron is in the state of discharge (black). At 20 s, the external stimulus is removed (turquoise), and the neuron cannot continue to fire and return (red) to the resting state (green). (b) At first, the neuron is in the state of discharge (black), At 20.01 s, the external stimulation is removed (turquoise), and the neuron continues to discharge and “delay” appeared (red). After that, “delay” ends (magenta) and the neuron returns (blue) to the resting state (green).

The values of td as ts varies from 0 s to 60 s are plotted in Fig. 4 for further studying the effect of ts on td. According to Fig. 4, four values of td can be clearly observed in the interval [20, 60] s for ts, which indicates that the firing of the neuron will only take on one of four patterns, with td values of 0 s, 2.8 s, 4.1 s or 5.1 s. Each of these four values is very dense, denoting frequent switching of td with changes in ts. The inserted graph in Fig. 4 reveals that these four values can be observed in an even narrower interval. It is suggested that this phenomenon has a fractal characteristic because of self-similarity. Therefore, the Hausdorff fractal dimension () can be used for quantification. Differing from the general calculation method, l is defined as the number of points chosen in Fig. 4 and N(l) is defined as the sum of the absolute values of the differences between the points chosen from l one-by-one.[47] The result of Hausdorff dimension calculation is shown in Fig. 5, which illustrates that this phenomenon does have fractal properties for its Hausdorff dimension is 0.9994. Therefore, when ts changes, even extremely subtly, td will be greatly altered.

Fig. 4. Delay (td) versus stimulus duration (ts) in the interval [0, 60] s of ts; there are 600001 points plotted in this picture. The inserted graph shows a detailed view of the interval [30, 31] s for ts.
Fig. 5. Hausdorff dimension calculation for the delay (td) versus stimulus duration (ts) plot as shown in Fig. 4, data in interval [20, 60] s are chosen for its clear and separate values.

In fact, there are a few scattered points with a value of approximately td = 1.8 s in the interval [20, 60] s in Fig. 4, which suggests that as the calculation precision improves, more points may appear and present as a horizontal line. In the ts interval [0, 20] s, more values can be observed than in the interval [20, 60] s; the distribution of these values may follow some kinds of rules. It may be that more regular values can be observed when the accuracy of the calculation is further improved.

3.2. The sensitive dependence of “delay” on the intensity of the stimulus

In addition to stimulus duration, the effect of stimulus intensity on neuronal firing and td is also studied. For these calculations, the stimulus duration is fixed at ts = 20 s.

Several neuron firing patterns are plotted in Fig. 6 under different stimulus intensities (Iext), the meaning of the two arrows is the same as that in Fig. 1. Figure 6(a) shows the neuron cannot begin to discharge without external stimulation, while the neurons cannot discharge normally although the external stimulation is enhanced (Figs. 6(b) and 6(d)). The values of td as a function of Iext from 0 μA/cm2 to 40 μA/cm2 are shown in Fig. 7.

Fig. 6. The sampled time series of neuronal firing patterns with ts = 20 s and different stimulus intensities. Neurons in panels (a), (b), and (d) cannot discharge well before the external stimulus ceases. The initial values are the same as those in Fig. 1.
Fig. 7. The delay (td) versus stimulus intensity (Iext) in the interval [0, 40] μA/cm2 of ts, the negative td means the stimulus intensity is not strong enough to activate that the neuron to fire 20 s at least; there are 40001 points in this picture.

The result in Fig. 7 does not truly indicate that td = −20 s in the Iext interval [0, 7.5] μA/cm2; according to our definition, a negative td is meaningless. As shown in Fig. 6(a), the neuron cannot fire continuously. Figures 6(b)6(d) show the firing patterns of the neuron in the interval [7.5, 15] μA/cm2, indicating different neuronal discharging modes that switch frequently. Some negative values of td reveal transient neuronal firing. In summary, the neuron either fails to fire or experiences a short period of time, less than 20 s, from the firing state to resting state when td is negative. The last interval, [15, 40] μA/cm2, shows results similar to those in Fig. 4. Moreover, there are primarily six values in the interval [7.5, 15] μA/cm2 and two values in the interval [15, 40] μA/cm2, which indicates the distribution of td is obviously related to the selection of interval. The results of the fractal dimension calculations in different intervals are shown in Fig. 8.

Fig. 8. Hausdorff dimension calculations for the delay (td) versus stimulus intensity (Iext) plot in the two intervals in Fig. 7, Data in two intervals are obtained separately, while the calculation method is the same as that in Fig. 4.

From Fig. 8, it can be seen that this phenomenon does have fractal characteristics. That is, Iext does cause severe switching in td, and the frequency of switching in different regions is different due to the fractal dimensions are 0.9837 for interval [7.5, 15] μA/cm2 and 0.966 for interval [15, 40] μA/cm2.

In the above two sections, the effects of ts and Iext towards td are studied respectively. However, both factors should be considered simultaneously for confirming the effect of this coupling system. Figure 9 shows that td is jointly influenced by these two factors; the gray and black areas indicate values of td outside the range shown on the right-hand side. Special intervals are chosen for clearer observation: ts is selected from 20 s to 60 s, and Iext is confined to [20, 40] μA/cm2. In the upper left corner, the value of td is outside the range shown on the right-hand side of the image. The largest value of td in this area is 17.295 s, and the smallest value is 0 s, with td changing more severely in this area. Moreover, changes along the horizontal axis are sharper than changes along the vertical axis, and the black-and-gray vertical bars also reflect this phenomenon. However, another interesting phenomenon should be noted. During the changes in the intensity and duration of the external stimulus, the value of td that appears most often is 5.1 s. Further research could clarify this feature of this coupling system. From Fig. 9, td shows severe switching, as a function of ts or Iext; therefore, high sensitivity is shown in both dimensions.

Fig. 9. The value of “delay” influenced by two factors: duration and intensity of stimulus. Values in gray and black colors are beyond the range shown on the right.

To investigate how the effect of electromagnetic induction regulates the neuron–astrocyte coupling system, the feedback strength k1 is varied to characterize changes in neuronal discharging modes and in td.

Figure 10 shows several distributions of td for different values of k1. Apparently, td increases with decreasing k1; as k1 approaches zero, the value of td increases further. Additionally, the values of td in the interval [60, 180] s of ts in Fig. 10(a) and in the interval [20, 60] s in the other three graphs of Fig. 10 are chosen because of its clearer distribution.

Fig. 10. The distributions of “delay” as a function of stimulus duration for different values of k1.

The value of td as k1 varies from 0 to 0.015 is plotted in Fig. 11(a); for this plot, the values of td are stable as described above. The exponential form of the plot implies a potential critical slowing down of the feedback strength. The critical value k1c = 0.004 is the most appropriate for matching the critical form: in all the values from 0 to 0.005; the linear fit of this equation in logarithmic coordinates is shown in Fig. 11(b). The successful fitting also proves the existence of the critical slowing down of k1.

Fig. 11. (a) The value of the stabilized “delay” (td) versus feedback strength k1; a potential critical slowing down of the feedback strength is implied in the exponential form of the plot. (b) Linear fitting of the log transformation of td and k1 according to the critical form (k1c = 0.004).

On the other hand, td will disappear when k1 is less than the critical value k1c = 0.004; it is certainly that there is no td in the absence of electromagnetic induction (k1 = 0). In other words, td appears precisely because of the existence of electromagnetic induction.

3.3. Hysteresis loop caused by changing the external stimulus intensity

In this section, the duration of the stimulus (ts) is set to be infinity, and the intensity of the stimulus (Iext) is varied. There are several neuronal discharge situations with astrocytic calcium ion concentrations under external stimulus intensities in the interval [7, 12] μA/cm2, as shown in Fig. 12. As the intensity of the external stimulus changes, the firing modes of the neuron and the regulation by the astrocyte seem to change irregularly. The sudden change of Iext from 8 μA/cm2 (Fig. 12(b)) to 9 μA/cm2 (Fig. 12(c)) leads to a big change of neuronal firing state, from resting to spiking, while it cannot fire normally at 10 μA/cm2 of the Iext (Fig. 12(d)). However, once the Iext reaches 11 μA/cm2 (Fig. 12(e)), the neuron can discharge continuously and the concentration of the astrocytic calcium ion is sustaining higher again.

Fig. 12. The firing patterns of the neuron and Ca2+ concentrations in the astrocyte with different stimulus intensities whose duration is extended to be infinite. The black curves represent neuronal membrane potentials, whereas the red curves represent calcium concentrations in the astrocytes. The initial values are the same as those in Fig. 1. (a) Iext = 7 μA/cm2, (b) Iext = 8 μA/cm2, (c) Iext = 9 μA/cm2, (d) Iext = 10 μA/cm2, (e) Iext = 11 μA/cm2, (f) Iext = 12 μA/cm2.

The concentration of Ca2+ is at a stable and high level when the neuron spikes continually, while it is at a low level when the neuron rests. Therefore, the concentration of Ca2+ in the astrocyte and the pattern of neuronal firing change simultaneously. The amplitude of the stable oscillations of calcium waves in the astrocyte is chosen as the ordinate because it is clearer when only one parameter is used as an indicator of the neuron firing pattern. As a result, the concentration of Ca2+ corresponding to Iext in the interval [0, 14] μA/cm2 is shown in Fig. 13 for studying its effect on the firing patterns of neurons; a stable Ca2+ concentration is set to the longitudinal coordinate, and the adjacent points are connected to indicate the severe change in Ca2+ concentration.

Fig. 13. Ca2+ concentration in a stable state versus stimulus intensity (Iext). The points are not continuous in the interval [8, 12] μA/cm2 and its vicinity, the adjacent parts are connected with lines for showing the sharp change in Ca2+ concentration.

From Fig. 13, the concentration of Ca2+ sharply switches with a change in Iext. Lower Ca2+ concentrations correspond to Figs. 12(a), 12(b), and 12(d), which is a stable focus in phase space. Higher concentrations are shown in Figs. 12(c), 12(e), and 12(f). The trajectories of the system in phase space is shown in Fig. 14 for more intuitive presentation. The magenta point and red circle represent the stable focus and the second attractor, respectively. There is a closed circle shown in Fig. 15, which is the Poincaré surface (Ca2+ = 315 nM) of the attractor for higher Ca2+ concentrations. The closed circle indicates that the second attractor of this system in phase space is an invariant torus attractor when higher Ca2+ concentrations are taken. That is to say, the whole system is in quasi-periodic motion. Although the initial values (blue–green) are the same, they will follow different trajectories (black or blue) to reach different attractors (magenta or red) under different external stimulus intensities. This is why the concentration of Ca2+ oscillates violently in the interval as depicted in Fig. 13.

Fig. 14. Two stable attractors and the trajectories leading to them; the initial value (blue–green) will follow different trajectories (black or blue) to reach different attractors (magenta or red) under different stimulus intensities.
Fig. 15. Poincaré surface of the invariant torus attractor for higher Ca2+ concentrations ([Ca2+] = 315 nM).

In fact, an increase in the external stimulus intensity causes the phase diagram of this multidimensional system to change in phase space. The initial value is switched between the attraction domain of the stable focus, whose time history is shown in Fig. 12(b), and the attraction domain of the invariant torus attractor, whose time history is shown in Fig. 12(c). The two attraction domains change as Iext changes.

In other words, when a value in the interval [7.5, 12.5] μA/cm2 is chosen for the external stimulus intensity, there will be at least two stable attractors in phase space; the neuronal discharging pattern and the oscillation of the concentration of Ca2+ have a bistable state between a quiescent state and quasi-periodic motion; the initial value may be taken near the boundaries of these two attractors, and the trajectory may be towards the focus or the torus attractor. Neuronal electric activity may consequently either be at rest or oscillate due to different initial states.

The initial values are the same in calculating the stable Ca2+ concentration under different Iext in Fig. 13. In order to find the accurate interval of bistability, the initial values at the next stimulus intensity is the values when the Ca2+ concentration converges under the present stimulus intensity, regardless of whether the intensity increases or decreases. Apparently, the initial values obtained above is different from the one shown in Fig. 13.

The entire bistable interval for Ca2+ concentrations is shown in Fig. 16, where a hysteresis loop is observed for an Iext from 3.51 μA/cm2 to 12.79 μA/cm2, which also illustrates the system’s sensitive selection of initial values.

Fig. 16. The bifurcation diagram of Ca2+ concentration versus external stimulus intensity (Iext); an intact bistable interval is shown over [3.51, 12.79] μA/cm2.

The effect of electromagnetic induction is considered as well; k1 is chosen as the variable once again to determine the dynamic effects of magnetic induction on astrocytic calcium ion concentrations.

Several bifurcation diagrams are shown in Fig. 17 for the bistable interval of Ca2+ concentration. According to Fig. 17, the lower critical external stimulus intensity (Il) and the higher critical external stimulus intensity (Ih) increase with increasing k1, while the difference between Il and Ih decreases. Moreover, the bistability disappears and the amplitude of the quasi-periodic oscillation becomes larger as k1 varies from 0.015 to 0.02.

Fig. 17. The bifurcation diagram of Ca2+ concentration with respect to external stimulus (Iext) under different k1, Il, and Ih mean lower and higher critical external stimulus of two catastrophe points.

To investigate how k1 specifically affects bistability, more values need to be calculated in the interval [0, 0.02]. Figure 18 shows the critical Iext as a function of k1. As the electromagnetic conduction increases from zero (i.e., k1 = 0), Il and Ih increase while the bistable interval is reduced. The gap between Il and Ih vanishes when k1 = 0.016, bistability disappears and only one catastrophe point is left. The magenta point is the critical value for the disappearance of bistability. Moreover, electromagnetic induction only acts on the whole system when it is strong, while it has little effect at a relatively lower level.

Fig. 18. The critical stimulus intensity (Iext) changes with the varying of k1, the magenta point indicates the moment bistability disappeared.
4. Conclusion and discussion

In this paper, we study the dynamic behavior of an improved neuron–astrocyte coupling model which takes into account the effect of electromagnetic induction.

The high sensitivity of the “delay”, the delayed time from the firing state to the resting state after the stimulus stopped, can be observed by changing the duration or intensity of the external stimulus. Although the control parameters are different, the effects of varying one parameter while keeping the other fixed are similar between the two parameters. This high sensitivity can even be shown in double-parameter space when both the duration and intensity are considered. The self-similarity implied in Fig. 4 indicates that several values can be taken even in very limited intervals, which can characterize the sensitivity. The fractal dimension can quantify the high sensitivity because self-similarity is also a characteristic of fractals. The results obtained from varying the stimulus duration or intensity also display the obvious fractal phenomenon. Moreover, different fractal characteristics can be observed over different intervals of the external stimulus intensity as the external stimulus intensity changes and are confirmed by the different calculated fractal dimensions. In addition, a bistable state about a quiescent state and quasi-periodic motion is observed for a certain range of stimulus intensities when the duration of stimulus is extended to infinity. The neuron could be in the resting state or the periodically oscillating state according to the initial values. Finally, the influence of electromagnetic induction on the coupled system was investigated. As the feedback strength (k1) of the electromagnetic induction on the neuronal membrane potential decreases to zero, the “delay” increases exponentially and is infinite in the critical value of k1 (k1 = 0.004) while it will disappear when k1 is decreasing further especially in the absence of electromagnetic induction (k1 = 0). That is, only electromagnetic induction contributes to the generation of “delay”. On the other hand, a strong electromagnetic induction will affect the bistability of the whole coupling system, but there is no evident effect of weak induction.

There are few researches on considering the magnetic effect of astrocytes, whether from the perspective of experiment or numerical study. However, magnetic effect or magnetic stimulation has a strong effect on neurocytes and neurological diseases. Repetitive transcranial magnetic stimulation (rTMS), an effective method in the treatment of depression for its anti-depressant effect associated with changes to the endocannabinoid system (ECS) and increase of the substance expression in astrocytes, has attracted attentions in recent years.[48,49] However, the magnetic stimulation put forward above is different from the magnetic effect caused by electric field induction in the present paper, which introduces the effect of magnetic induction to neurons coupled with astrocytes at the same time.

For some time, it had been believed that electromagnetic induction contributed to the memory effect of the nervous system,[34,50] which is typically realized by introducing a time delay term into the model. Obviously, the time delay is different from the “delay” defined here, which is the delayed duration back to resting state while the external stimulation ceased in the neuron–astrocyte coupling model that takes into account electromagnetic induction. Moreover, the effects of sustained external stimulation on neurons have been considered in Refs. [41,42], while “delay” only occurs when there is limited stimulation duration according to the present paper. A magnetic matter in the brain, which was proved to be the magnetic protein called MagR, was successful predicted due to the inductance effect by energy coding model,[51] it provides a perspective for exploring further. The model defined here may create a novel viewpoint for studying the effects of neuronal systems.

In the clinic, numerous methods involving electric stimulation are used to treat nervous system diseases.[52] Vagus nerve stimulation therapy has been shown to have an effect in the treatment of epilepsy,[53] and electrical stimulation of the cervical or thoracic spinal cord can significantly help patients with dystonia and spastic torticollis to control motor function.[54] However, depending on the type of disease and patient tolerance, subtle changes in stimulus intensity and duration could have a significant impact on the therapeutic effect. The sensitivity illustrated in this paper may provide a new perspective for explaining this phenomenon.

Reference
[1] Hodgkin A L Huxley A F 1952 J. Physiol. 117 500
[2] Morris C Lecar H 1981 Biophys. J. 35 193
[3] Fitzhugh R 1960 J. Gen. Physiol. 43 867
[4] Hindmarsh J L Rose R M 1982 Nature 296 162
[5] Zhao Z G Gu H G 2017 Sci. Rep. 7 1
[6] Bo L Liu S Q Liu X L et al. 2016 Int. J. Bifur. Chaos 26 1650090
[7] Gu H G Zhao Z G Bing J et al. 2015 Plos One 10 e0121028
[8] Wang H X Wang Q Y Zheng Y H 2014 Sci. China: Technol. Sci. 57 872
[9] Zhao Z G Bing J Gu H G 2016 Nonlinear Dyn. 86 1549
[10] Liu X L Liu S Q 2012 Nonlinear Dyn. 67 847
[11] Gu H G Pan B B Li Y Y 2015 Nonlinear Dyn. 82 1191
[12] Ma J Wu F Q Alsaedi A et al. 2018 Nonlinear Dyn. 93 2057
[13] Gu H G Zhao Z G 2015 Plos One 10 e0138593-
[14] Ma J Zhang G Hayat T et al. 2019 Nonlinear Dyn. 95 2019
[15] Wang R B Wang Z Y Zhu Z Y 2018 Nonlinear Dyn. 92 973
[16] Wang R B Ichiro T Zhang Z K 2015 Int. J. Neural Syst. 25 1450037
[17] Zhu F Y Wang R B Pan X C 2019 Cognitive Neurodynamics 13 75
[18] Wang Z Y Wang R B 2014 Frontiers in Computational Neuroscience 8 14
[19] Zhu Z Y Wang R B Zhu F Y 2018 Front. Neurosci. 12 122
[20] Wang Y H Wang R B Zhu Y T 2017 Cognitive Neurodynamics 11 99
[21] Wang Y H Xu X Y Wang R B 2018 Front. Neurosci. 12 264
[22] Porter J T Mccarthy K D 1996 J. Neurosci. 16 5073
[23] Parpura V Haydon P G 2000 Proc. Natl. Acad. Sci. USA 97 8629
[24] Arcuino G Lin J H Takano T et al. 2002 Proc. Natl. Acad. Sci. USA 99 9840
[25] Postnov D E Ryazanova L S Brazhe N A et al. 2008 J. Biol. Phys. 34 441
[26] Li Y X Rinzel J 1994 J. Theor. Biol. 166 461
[27] Nadkarni S Jung P 2003 Phys. Rev. Lett. 91 268101
[28] Li J J Xie Y Yu Y G et al. 2017 Sci. China: Technol. Sci. 60 43
[29] Erkan Y Sarac Z Yilmaz E 2019 Nonlinear Dyn. 95 3411
[30] Kanakov O Gordleeva S Ermolaeva A et al. 2019 Phys. Rev. 99 012418
[31] John F B Matthew J T Jennifer M S et al. 2016 Proc. Natl. Acad. Sci. 113 14133
[32] Alex M A 2002 Kybernetes Emerald Group Publishing Limited 140 142 10.1108/03684920210413818
[33] Ma J Tang J 2017 Nonlinear Dyn. 89 1569
[34] Lv M Wang C N Ren G D et al. 2016 Nonlinear Dyn. 85 1479
[35] Zhang X H Liu S Q 2018 Chin. Phys. 27 040501
[36] Chua L 1971 IEEE Trans. Circuit Theory 18 507
[37] Wang R B Zhang Z K Jiao X F 2006 Appl. Phys. Lett. 89 1102
[38] Wang Z Y Wang R B Fang R Y 2015 Cognitive Neurodynamics 9 129
[39] Wang R B Zhu Y T 2016 Cognitive Neurodynamics 10 1
[40] Feng P H Wu Y Zhang J Z 2017 Frontiers in Computational Neuroscience 11 94
[41] Liu Y Ma J Xu Y et al. 2019 Int. J. Bifur. Chaos 29 1950156
[42] Wu F Q Gu H G Li Y Y 2019 Commun. Nonlinear Sci. Numer. Simul. 79
[43] Muramatsu S Toda M Nishikawa J et al. 2019 Brain Res. 1721
[44] Katta S Sanzeni A Das A et al. 2019 J. Gen. Physiol. 151 1213
[45] Wang S S Alousi A A Thompson S H 1995 J. Gen. Physiol. 105 149
[46] Ma J Tang J 2015 Sci. China: Technol. Sci. 58 2038
[47] Budyansky M Uleysky M Prants S 2002 Physica D-Nonlinear Phenomena 195 369
[48] Xue S S Xue F Ma Q R et al. 2019 Pharmacol., Biochem. Behav. 184
[49] Zorzo C Higarza S G Mendez M et al. 2019 Brain Res. Bull. 150 13
[50] Wang R Feng P H Fan Y C et al. 2019 Int. J. Bifur. Chaos 29 1950005
[51] Wang Y Y Wang R B 2018 Nonlinear Dyn. 91 319
[52] Kerrigan J F Litt B Fisher R S et al. 2010 Epilepsia 45 346
[53] Steven C S Clifford B S 1998 Epilepsia 39 677
[54] Dooley D M Nisonson I 1981 Appl. Neurophysiol. 44 71